
library(plm);library(countrycode);library(stargazer)


data<-read.csv(file="data_final.csv")


dist_matrix<-read.csv("dist_matrix2.csv")


dist_matrix$X<-countrycode(dist_matrix$X,"country.name","cown")

names(dist_matrix)[2:ncol(dist_matrix)]<-dist_matrix$X


data<-data[data$year > 1940,]



dd<-c(.4234/2,.4234,.4234*2)


reg_out2<-reg_out<-list()


for(k in 1:length(dd)){
  
  
  out<-list()
  
  for(i in 1:nrow(dist_matrix)){
    
    out2<-list()
    
    for(j in min(data$year):max(data$year)){
      
      sub<-data[data$ccode==dist_matrix$X[i] & data$year==j, ]
      
      if(nrow(sub) > 0){
        
        sub_year<-data[data$year==j,]
        
        sub_y<-data.frame(sub_year$y,sub_year$D,sub_year$ccode)
        
        names(sub_y)<-c("y","D","ccode")
        
        
        sub_y<-sub_y[order(sub_y$ccode),]
        
        
        dist1<-dist_matrix[dist_matrix$X==sub$ccode,]
        
        dist1<-dist1[1,]
        
        
        dist2<-dist1[names(dist1)%in%sub_y$ccode]
        
        dist2<-dist2[names(dist2) %in% sub_y$ccode]
        
        sub_y<-sub_y[sub_y$ccode %in% names(dist2), ]
        
        
        sub_y$dist<-t(as.vector(dist2))
        sub_y$dist<-sub_y$dist/1000000
        
        sub_y<-sub_y[sub_y$dist!=0,]
        
        
        
        sub$Y_dist_D<-sum( sub_y[!is.na(sub_y$y)& sub_y$D==1,]$y * exp(-sub_y[!is.na(sub_y$y)& sub_y$D==1,]$dist*dd[k]) )  / sum(exp(-sub_y[!is.na(sub_y$y)& sub_y$D==1,]$dist*dd[k] ))
        
        
        sub$Y_dist_A<-sum( sub_y[!is.na(sub_y$y)& sub_y$D==0,]$y * exp(-sub_y[!is.na(sub_y$y)& sub_y$D==0,]$dist*dd[k] ) )  / sum(exp(-sub_y[!is.na(sub_y$y)& sub_y$D==0,]$dist*dd[k] ))
        
        
        sub$D_dist<-sum( sub_y[!is.na(sub_y$y),]$D * exp(-sub_y[!is.na(sub_y$y),]$dist*dd[k] ) )  / sum(exp(-sub_y[!is.na(sub_y$y),]$dist*dd[k] ))
        
        
        
        
        
      }
      
      
      if(nrow(sub)==0){
        
        sub<-rep(NA,14)
        
        names(sub)<-c(names(data),"Y_dist_D","Y_dist_A","D_dist")
        
      }
      
      
      
      
      out2[[j-min(data$year)+1]]<-sub
      
      
    }
    
    out[[i]]<-do.call(rbind,out2)
    
    
    
    
  }
  
  
  ldata<-do.call(rbind,out)
  
  ldata<-ldata[!is.na(ldata$ccode),]
  
  ldata<-ldata[ldata$ccode != 345,]
  
  ldata<-ldata[ldata$year > 1947 & ldata$year <=2000,]
  
  
  pdata<-pdata.frame(ldata,index=c("ccode","year"))
  
  pdata$Y_diff<-pdata$Y_dist_D-pdata$Y_dist_A
  
  pdata$l3yda<-(lag(pdata$Y_dist_A,1)+lag(pdata$Y_dist_A,2)+lag(pdata$Y_dist_A,3))/3 
  pdata$l3ydd<-(lag(pdata$Y_dist_D,1)+lag(pdata$Y_dist_D,2)+lag(pdata$Y_dist_D,3))/3 
  
  
  pdata$l3ydiff<-(lag(pdata$Y_diff,1)+lag(pdata$Y_diff,2)+lag(pdata$Y_diff,3))/3 
  
  
  pdata$l3ddf<-(lag(pdata$D_dist,1)+lag(pdata$D_dist,2)+lag(pdata$D_dist,3))/3
  
  pdata$ld<-lag(pdata$D,1)
  
  pdata$ly<-lag(pdata$y,1)
  
  
  pdata$l3y<-(lag(pdata$y,1)+lag(pdata$y,2)+lag(pdata$y,3))/3
  
  
  #pdata$fail<-ifelse(pdata$TinP ==0,1,0)
  
  #pdata$fail<-1-pdata$I
  
  
  
  pdata$fail<-pdata$I
  
  
  pout<-plm(D~ld+l3ddf+l3ydiff,model="within",data=pdata)
  
  pout1<-plm(D~ld+l3ddf+l3ydiff+log(lag(Y,1)),model="within",data=pdata)
  
  pout2<-plm(D~ld+l3ddf+l3ydiff+TinP+log(lag(Y,1)),model="within",data=pdata)
  
  pout3<-plm(D~ld+l3ddf+l3ydiff+TinP+log(lag(Y,1))+lag(Trade,1),model="within",data=pdata)
  
  
  pout4<-plm(fail~D*y,model="within",data=pdata)
  
  
  pout4a<-plm(fail~D*y | l3ddf+ l3yda+l3ydd,model="within",data=pdata)
  
  pout5<-plm(fail~D*y+log(lag(Y,1)),model="within",data=pdata)
  
  
  pout5a<-plm(fail~D*y+log(lag(Y,1)) | l3ddf+ l3yda+l3ydd+log(lag(Y,1)),model="within",data=pdata)
  
  
  
  pout6<-plm(fail~D*y+TinP+log(lag(Y,1)),model="within",data=pdata)
  
  pout6a<-plm(fail~D*y+log(lag(Y,1)) +TinP | l3ddf+ l3yda+l3ydd+log(lag(Y,1))+TinP,model="within",data=pdata)
  
  
  pout7<-plm(fail~D*y+TinP+log(lag(Y,1))+lag(Trade,1),model="within",data=pdata)
  
  pout7a<-plm(fail~D*y+log(lag(Y,1)) +TinP +lag(Trade,1)| l3ddf+ l3yda+l3ydd+log(lag(Y,1))+TinP+lag(Trade,1),model="within",data=pdata)
  
  
  
  
  reg_out[[k]]<-list(pout,pout1,pout2,pout3)
  
  
  reg_out2[[k]]<-list(pout4a,pout5a,pout6a,pout7a)
  
  
  reg_out3<-list(pout4,pout5,pout6,pout7)
  
}





library("stargazer")


### Table A2


stargazer(list(reg_out[[1]],reg_out[[2]],reg_out[[3]]), keep= c("ld","l3ydiff","l3ddf"))



### Table A3


stargazer(list(reg_out3,reg_out2[[2]]),keep=c("D","y","D:y"))






##############  Figure A 1  ####################################3


data$D_year<-paste(data$year,data$D,sep="-")


d_diff<-tapply(data$y,data$D_year,mean,na.rm=T)



d_diff<-data.frame(names(d_diff),d_diff)


names(d_diff)<-c("D_year","y")


d_diff$year<-do.call(rbind,strsplit(as.character(d_diff$D_year),"-"))[,1]
d_diff$D<-do.call(rbind,strsplit(as.character(d_diff$D_year),"-"))[,2]








d_diff<-d_diff[as.numeric(as.character(d_diff$year)) > 1949,]






pd<-paste(wd, "/Appendix_A2/A_output/",sep="")


pdf(file= paste(pd,"diff_growth_rates.pdf",sep=""), 12, 6)


par(mar=c(3,5,2,1))

plot(d_diff[d_diff$D==1,"year"],d_diff[d_diff$D==1,"y"]-d_diff[d_diff$D==0,"y"],type="l",lwd=2,col="black",ylim=c(-.015,.035),xlim=c(1950,2000),axes=F,xlab="",ylab=expression(bar(y)^{D==1}-bar(y)^{D==0}) )
grid()
axis(1);axis(2)
abline(h=0,lty=2,col="red")

mtext("Difference in Observed Growth Rates",cex=1.5,font=2)

par(new=T)


D_diff<-tapply(data$D,data$year,mean,na.rm=T)

D_diff<-data.frame(names(D_diff),D_diff)

names(D_diff)<-c("year","D")

D_diff$year<-as.numeric(as.character(D_diff$year))

D_diff<-D_diff[D_diff$year > 1949,]

D_diff$diff<-d_diff[d_diff$D==1,"y"]-d_diff[d_diff$D==0,"y"]




dev.off()


